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General relativistic simulations of binary black hole-neutron stars: 
Precursor electromagnetic signals 

Vasileios Paschalidis, Zachariah B. Etienne, and Stuart L. Shapirqj 
Department of Physics, University of Illinois at Urbana- Champaign, Urbana, IL 61801 

We perform the first general relativistic force-free simulations of neutron star magnetospheres in 
orbit about spinning and non-spinning black holes. We find promising precursor electromagnetic 
emission: typical Poynting luminosities at, e.g., an orbital separation of r = 6.6i?NS are Lem ~ 
6 x lO 42 ( J B N s,p/lO 13 G) 2 (MNs/1.4M ) 2 erg/s. The Poynting flux peaks within a broad beam of 
~ 40° in the azimuthal direction and within ~ 60° from the orbital plane, establishing a possible 
lighthouse effect. Our calculations, though preliminary, preview more detailed simulations of these 
systems that we plan to perform in the future. 

PACS numbers: 04.25. D-, 04.25. dk,04.30.-w,52. 35. Hr 



Black hole-neutron star (BHNS) binaries are among 
the most promising sources for the simultaneous detec- 
tion of gravitational wave (GW) and electromagnetic 
(EM) signals in the era of multimessenger astronomy. For 
example, Advanced LIGO is expected to detect between 
1-100 BHNS GW signals each year 0-Q. In addition, 
BHNS mergers may provide the central engine that pow- 
ers a short-hard gamma-ray burst (sGRB). GW signals 
from the inspiral and merger of BHNSs have been com- 
puted recently in full general relativity (GR) [5l4ll|, and 
the first parametric study of magnetized BHNS mergers 
in full GR has been been carried out in [12|, [l3[ , where it 
was shown that under appropriate conditions BHNSs can 
launch collimated jets - necessary ingredients for many 
sGRB models. 

Apart from post-merger EM signals, detectable EM 
pre-merger signals may arise during a BHNS inspiral. 
Detecting such precursor signals, combined with GW 
observations, will yield a wealth of information about 
BHNS binaries. EM signals will help localize the source 
on the sky, resulting in improved parameter estimation 
from GWs by eliminating degeneracies resulting form im- 
precise localization of the source [14]. 

Neutron stars are believed to be endowed with dipole 
magnetic fields. Thus, toward the end of a BHNS inspi- 
ral, strong magnetic fields will thread and sweep the BH 
horizon, possibly establishing a unipolar inductor (UI) 
DC circuit [15| that extracts energy from the system. 
This exciting new possibility has been suggested recently 
[16j as a potential mechanism for powering precursor EM 
signals from BHNSs. Follow-up analytical approxima- 
tions in the high-mass ratio limit have been performed 
16l~ll8j to estimate the total EM output. But, as these 
UIs operate in strongly curved, dynamical spacetimes, 
numerical relativity simulations are necessary to reliably 
determine the amount of EM output, especially in the 



regime of comparable-mass binaries where previous ap- 
proximations do not apply. While UIs are not unique to 
BHNS binaries (they can operate in other types of com- 
pact binaries, such as NSNSs 17|, [l9|, [20| ) , BHNSs may 
be optimal systems for this mechanism because the az- 
imuthal twist of the magnetic flux tubes is less than unity 
for a BH resistor 
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In this Letter we simulate NS magnetospheres in or- 
bit about spinning and nonspinning BHs prior to merger 
via general relativistic, force- free (GRFF) simulations. 
We calculate the Poynting luminosity and characterize 
its angular dependence. In addition, we treat another 
EM emission mechanism: magnetic dipole (MD) radi- 
ation from the accelerating NS. MD radiation has been 
considered in the context of EM emission affecting the in- 
spiral and GW signal [22|, but not as a potential source 
for strong precursor EM signals. Here we show that the 
MD Poynting luminosity is significant, and may be the 
dominant source of EM output in cases where UI does 
not operate either because of corotation or large values 
of azimuthal twist of the magnetic flux tubes. We use 
geometrized units where c = 1 = G, unless otherwise 
stated. 

Solving the GRFF equations generally involves evolv- 
ing the electric (E) and magnetic (B) fields under the 
force-free constraints E • B = and E 2 < B 2 [H, Q. 
The force-free regime represents the limit of ideal MHD 
when the magnetic fields dominate the plasma dynam- 
ics [23|, [25|- In this regime one can choose the B-field 
and the Poynting vector S as the dynamical variables 
and cast their evolution equations in conservation form 
[26, 27]. The force- free constraints then become S-B = 
and S 2 < B A as we show in [27j . One of the greatest ad- 
vantages of this alternative GRFF formulation is that an 
ideal GRMHD code can be easily modified to also solve 
the GRFF equations [26]. The GRFF formulation we 
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adopt here is the same as that in [26| , except that at ev- 
ery timestep, in addition to S 2 < £? 4 , we also enforce the 
algebraic constraint S • B = 0, which was ignored in [26|. 
This formulation is embedded in the fully GRMHD in- 
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FIG. 1. Initial magnetic field in the nonspinning a* = case (upper left panel). Relaxed magnetic field at t « 1.5 orbits: 
a* = (lower left panel), a* = —0.5 (upper right panel), and a* = 0.75 (lower right panel). The black sphere represents the 
BH horizon and the NS is shown in red. Both white and yellow lines are the magnetic fields lines. White lines distinguish field 
lines that intersect the BH horizon. 



frastructure we presented and tested in [28l-l30j . In addi- 
tion, to enforce the V • B = constraint on our adaptive- 
mesh-refinement grids, the magnetic induction equation 
is evolved via the vector potential formulation we intro- 
duced in [12|, l29|, [30| , coupled to the Generalized Lorenz 
(GL) gauge condition we devised [13|, |30|, [3l| . We set the 
damping parameter of GL gauge £ = 1.5/ At, where At 
is the timestep of the coarsest level. 

At large separations, the inspiral timescale is much 
longer than the orbital timescale. So to model the 
BHNS binary spacetime and the NS matter fields, we 
adopt quasiequilibrium solutions of the conformal-thin- 
sandwich (CTS) equations for companions at fixed or- 
bital separation [5|, |32|, [33| . The CTS approximation is 
excellent at the separations and BH spins we consider in 
this work. It also yields a binary spacetime with a heli- 
cal Killing vector, which enables us to perform the sim- 
ulations in the center-of-mass frame by simply rotating 
the spacetime and the fluid rest-mass density and four- 
velocity in the same fashion as in [34]. This implemen- 
tation is equivalent to assuming stationary equilibrium 
of the matter and the gravitational fields in the corotat- 
ing frame of the binary. The problem is thus reduced 
to evolving the EM fields (B and S) in the background 
matter fields and spacetime. 

Given that force- free electrodynamics is a limit of ideal 
MHD, the same ideal MHD evolution equations can be 
used to evolve both the NS interior and the force-free 



exterior EM fields, provided in the exterior a compatible 
force-free velocity is used [26] and the rest-mass density 
is set to zero. This guarantees a smooth transition from 
the ideal MHD interior to the force- free exterior, and 
the MHD variables on the NS surface effectively provide 
boundary conditions for the exterior force-free evolution. 
However, given that the chosen initial A- field is not a 
CTS solution, we evolve the induction equation [Eqs. (8), 
(9) in 12]] in the NS interior, using the known CTS fluid 
four- velocity. This sets the boundary condition on the 
NS surface for the force-free velocity and magnetic field 
in the exterior. For more details see [27j. Note that after 
tidal disruption, a GRFF treatment becomes inadequate 
and must be replaced by full GRMHD. 

In addition to our new GRFF evolution techniques, 
we have also added two new (but equivalent) diagnos- 
tics to monitor the outgoing EM luminosity: (i) the <j)2 
Newman-Penrose scalar [3514371], and (ii) the Poynting 
vector S = (E x B)/47r. To compute <j)2 we use the same 
null tetrad as in [38| , and the outgoing luminosity is [39] 



Lem = lim 
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r 2 S^dn. (1) 



The spacetime and NS initial data we use correspond 
to cases A, B, C given in Table I in [5]. The BH spin pa- 
rameters are a* = cl/Mh = —0.5, 0, 0.75, and the BH:NS 
mass ratio is q — 3. The CTS NS fluid is modeled as 
an equilibrium, irrotational, unmagnetized, T = 2 poly- 
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FIG. 2. Angular distribution of Poynting flux, normalized by its peak value on a sphere of radius 120M = 915(MNs/l-4M0)km. 
Left: spin -0.5, middle: spin 0, right: spin 0.75. The plots correspond to a time after about 2 orbits. 



trope. We seed the initial NS with a purely poloidal mag- 
netic field that approximately corresponds to that gen- 
erated by a current loop. The coordinate-basis toroidal 
component of this vector potential is 



Aa> = 



r b 2 o 



( r 2 +r 2)3/2 



1 



15r 2 (r 2 +^ 2 ) 
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where vq is the current loop radius, Iq the loop current, 



r 2 = (x-x NS ) 2 + (y-yNsr + z 
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2/ns) 2 , and £ns,2/ns are the initial coordinates of the NS 
center of mass. For roCr Eq. (J2j) gives rise to the stan- 
dard B-field from a current loop on the z-axis, and the 
characteristic 1/r 3 fall-off of a standard magnetic dipole 
on the z = plane. Choosing ro = i?Ns/3 in all our sim- 
ulations, where R^s is the NS polar radius, we find that 
the initial magnetic field scales as 1/r 3 outside the NS to 
a very good degree. Our simulations are independent of 
the magnetic field strength. If we set Iq = 0.0007, the 
initial NS polar magnetic field (as measured by a CTS 
normal observer) is 8.8 x 10 15 G. The initial magnetic field 
geometry is shown in the upper left panel of Fig. HJ 

For the a* = case we perform a resolution study: 
the low, medium and high resolutions cover the BH hori- 
zon radius by 19, 29 and 36 zones, respectively. We also 
perform the simulations of the a* ^ BH cases using 
two resolutions, corresponding to the medium and high 
resolutions of the a* = case. In all simulations we 
use 9 levels of refinement placing the outer boundary 
at 400M « 3O5O(M NS /1.4M )km, and impose reflection 
symmetry across the orbital plane. 

After a transient phase that lasts a little over 1 orbit, 
the magnetic field geometry settles into a quasistationary 
configuration, which we show in Fig. [TJ It is evident that 
in the spinning cases, partial winding of the magnetic 
field has taken place due to frame dragging, which is 
most prominent in the a* =0.75 case. 

In Fig. [2]we show the angular distribution of the Poynt- 
ing flux. In all cases, the Poynting flux peaks within a 
broad beam of ~ 40° in the azimuthal direction, and in 
the a* = and a* = 0.75 cases within ~ 60° from the 
orbital plane. This may establish a lighthouse effect as 



a characteristic EM signature of BHNS systems prior to 
merger, if the variation is not washed out by intervening 
matter. 

The time evolution of the computed luminosities is 
shown in Fig. [3l After a transient period caused by our 
choice of non-stationary initial magnetic fields, the lu- 
minosities settle to an approximately constant value as 
expected. We find that the time-averaged luminosities 
after the first 1.5 orbits at the adopted separation are 
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(3) 

where 5ns, p is the NS polar magnetic field strength mea- 
sured by a CTS normal observer, and Mns is the NS rest 
mass. As the magnetic field inside the NS does not feed 
back onto the matter evolution, the EM luminosity scales 
exactly as B 2 for any field strength. The characteristic 
frequency of this EM radiation is of order the orbital fre- 
quency ~ 2OO(Mns/1-4M0) _1 Hz at the adopted separa- 
tion, and hence smaller than typical interstellar-medium 
plasma frequencies ~ 9kHz. Thus, this radiation will be 
reprocessed before it reaches the observer. 

We now compare our results to the approximate UI 
formula. The Poynting luminosity of a BHNS UI in the 
large q limit is given by [16] 
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where rn is the horizon radius in units of the BH mass 
Mh, 5ns, p is the NS polar magnetic field as measured 
by zero-angular- momentum observers (ZAMOs) [40], i.e., 
normal observers in a Kerr spacetime in Boyer-Lindquist 
coordinates, and r is the binary separation. Here v re \ is 
the azimuthal velocity of magnetic field lines as measured 
by ZAMOs, for which the following relation was proposed 
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FIG. 3. Poynting luminosity as a function of time calculated 
on a sphere of radius 120M = 915(M N s/1.4M )km for all 3 
cases: a* = —0.5 (red) dashed line, a* = (black) solid line, 
a* = 0.75 (magenta) dotted line. The inset focuses on the 
last 1.7 orbits of evolution. Here B p ,i3 — Bns, p /10 13 G and 
T rb is the orbital period. 



[iflj: v re} = r(Q - Q N s) 



where Q is the orbital 



'ns; - 4^, 

angular frequency, and Q^s is the NS spin angular fre- 
quency. Given that our BHNS binaries are irrotational, 
we set Qns = 0. Using the binary parameters from our 
simulations and setting 5ns, p ~ ^ns,p i n Eq. (|1]), we 
find 



£ui,a*=0.75 

^UI,a*=0 

£ui,a*=-0.5 



0.12(L o . =0 .75>, 

0.5(L a . =0 ), 
0.7(L , = _o.5>. 



(5) 



Thus, the UI formula seems to predict well the 
overall magnitude of our computed luminosi- 
ties. However, in contrast to Eq. ([5]), which 
predicts that £ui,a*=-o.5/£ui,a*=o ~ 1.5 and 
^ui,a*=-o.5/^ui,a,=o.75 « 7.8, our calculations © 
show only a weak dependence of the Poynting luminosity 
on the BH spin. This is likely due in part to the spin 
dependence being added linearly in the proposed formula 
for ?; re i, and in part to the existence of magnetic dipole 
emission. 

In addition to the UI luminosity, another important 
EM radiation emission mechanism, that always operates, 
is that due the accelerating MD moment of the NS. This 
has been considered previously [22| in the context of EM 
interactions as they affect the binary inspiral, but not 
as a mechanism for strong precursor EM counterparts to 
GWs. The approximate MD luminosity is given by [22j 
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FIG. 4. Convergence of EM luminosity normalized by the 
maximum luminosity vs. time. The difference between high 
and medium resolutions is smaller than that between medium 
and low resolutions, indicating that our scheme is convergent. 

where we inserted parameters from our simulations. 
Lem,md is only ~ 20 times smaller than what we observe 
in our simulations, but is included in our calculated lumi- 
nosity. MD emission dominates when UI does not oper- 
ate either because of corotation (and a* = for BHNSs) 
or because it cannot be established due to the azimuthal 
twist of the magnetic flux tubes being too large, which 
may be the case for NSNS binaries [21] . 

The results of our resolution study in the a* = case 
are shown in Fig. HJ where it is demonstrated that our 
scheme is convergent and that the resulting luminosity in 
the two highest resolutions agree to within ~ 5%. Due to 
numerical resistivity, the EM energy in the NS interior 
is conserved after 3 orbits to within 10%, 11%, 7% in 
the a* = —0.5, a* = 0, a* = 0.75 cases, respectively. 
Thus, these errors should be taken as the approximate 
error bars of our calculations. 

In a future work we plan to extend our simulations to 
study the variation of the outgoing Poynting luminosity 
during the inspiral phase, and its dependence on different 
mass ratios. 
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